function [ImageBinary,coveragePercentage]=PlotSynergy_L145_limb_Apex(L1,L4,L5,epoch,FN_i,CameraAngleCapability,DU,VU,TU,h_apex)

options=odeset('RelTol',1e-13,'AbsTol',1e-13);
R_sun = 696340;
theta = asind(R_sun/DU);
alpha = acosd(R_sun/(R_sun+h_apex));
beta = acosd(R_sun/(R_sun+h_apex));
Theta_min = 90-alpha-theta;
Theta_max = 90+beta-theta;

%% On-disk
temp_time_month = month(L4.epoch);
ind = find(temp_time_month==epoch);
[row,col,~] = size(FN_i);
visibility_disk=zeros(row,col);
visibility_limb=zeros(row,col);


S_i = L4.S_i;
for ii = 1:row
    for jj = 1:col
        normvec = [FN_i(ii,jj,1);FN_i(ii,jj,2);FN_i(ii,jj,3)];
        for k = ind(1)
            posvec = S_i(1:3,k);
            angle_temp = acosd(dot(normvec,posvec)/(norm(normvec)*norm(posvec)));
            if angle_temp < CameraAngleCapability
                visibility_disk(ii,jj) = 1;
            end
            if and((Theta_min<angle_temp),(angle_temp<Theta_max))
                visibility_limb(ii,jj) = 2;
            end
            if round(angle_temp)==0
                L4Location=[ii,jj];
            end
        end
    end
end
S_i = L5.S_i;
for ii = 1:row
    for jj = 1:col
        normvec = [FN_i(ii,jj,1);FN_i(ii,jj,2);FN_i(ii,jj,3)];
        for k = ind(1)
            posvec = S_i(1:3,k);
            angle_temp = acosd(dot(normvec,posvec)/(norm(normvec)*norm(posvec)));
            if angle_temp < CameraAngleCapability
                visibility_disk(ii,jj) = 1;
            end
            if and((Theta_min<angle_temp),(angle_temp<Theta_max))
                visibility_limb(ii,jj) = 2;
            end
            if round(angle_temp)==0
                L5Location=[ii,jj];
            end
        end
    end
end
S_i = L1.S_i;
for ii = 1:row
    for jj = 1:col
        normvec = [FN_i(ii,jj,1);FN_i(ii,jj,2);FN_i(ii,jj,3)];
        for k = ind(1)
            posvec = S_i(1:3,k);
            angle_temp = acosd(dot(normvec,posvec)/(norm(normvec)*norm(posvec)));
            if angle_temp < CameraAngleCapability
                visibility_disk(ii,jj) = 1;
            end
            if and((Theta_min<angle_temp),(angle_temp<Theta_max))
                visibility_limb(ii,jj) = 2;
            end
            if round(angle_temp)==0
                earthLocation=[ii,jj];
            end
        end
    end
end

visibility=visibility_disk+visibility_limb;


ImageBinary=mat2gray(visibility);
Window = figure('OuterPosition',[100,100,900,500]);
hold on;
mesh(ImageBinary,'EdgeColor','interp','FaceColor','interp');
xlabel('\lambda [deg]')
ylabel('\phi [deg]')
titlestring = sprintf('L4 Month = %d',epoch);
% title(titlestring)
set(gca,'fontsize',15)
xticks(linspace(0,row,5));
xticklabels({'180W','90W','0','90E','180E'});
yticks(linspace(0,col,5))
yticklabels({'90S','45S','0','45N','90N'});
xlim([0,row]);
ylim([0,col])
colormap("turbo");

counter = 0;
for ii = 250-111:250+111
    for jj = 1:col
        if visibility(ii,jj) == 3
            counter=counter+1;
        end
    end
end

plot3(earthLocation(2),earthLocation(1),1000,'ro','MarkerSize',15,'MarkerFaceColor','r')
plot3(L4Location(2),L4Location(1),1000,'rs','MarkerSize',15,'MarkerFaceColor','r')
plot3(L5Location(2),L5Location(1),1000,'rd','MarkerSize',15,'MarkerFaceColor','r')

coveragePercentage = counter/(row*col*(80/180));
fprintf('Percent Coverage : %f percent \n\n',coveragePercentage*100)


figure;
hold on;
[xs,ys,zs] = ellipsoid(-0,0,0,696340,696340,696340,500);
figure;
s = surface(xs,ys,zs);
rotate(s, [1 0 0], 7.25);
rotate(s, [0 0 1], 73.67+0.014*(2023-1850));
InclinationRotationAngle = 7.25;
RAANRotationAngle = 73.67+0.014*(2023-1850);
PoleDirection = Rotation([0;0;1],RAANRotationAngle,'Degrees')*...
                Rotation([1;0;0],InclinationRotationAngle,'Degrees')*...
                [0;0;1];
rotate(s, PoleDirection, 0);
surface(xs,ys,zs,ImageBinary_L145,'EdgeColor','none','AlphaData',0.5,'FaceAlpha',0.6);
axis equal
end